Multifractal finite-size-scaling and universality at the Anderson transition 
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We describe a new multifractal finite size scaling (MFSS) procedure and its application to the 
Anderson localization-delocalization transition. MFSS permits the simultaneous estimation of the 
critical parameters and the multifractal exponents. Simulations of system sizes up to L^ — 120"^ 
and involving nearly 10^ independent wavefunctions have yielded unprecedented precision for the 
critical disorder Wc = 16.530(16.524, 16.536) and the critical exponent u = 1.590(1.579, 1.602). We 
find that the multifractal exponents A, exhibit a previously predicted symmetry relation and we 
confirm the non-parabolic nature of their spectrum. We explain in detail the MFSS procedure first 
introduced in our Letter [Phys. Rev. Lett. 105, 046403 (2010)] and, in addition, we show how to take 
account of correlations in the simulation data. The MFSS procedure is applicable to any continuous 
phase transition exhibiting multifractal fluctuations in the vicinity of the critical point. 

PACS numbers: 71.30.+h,72.15.Rn,05.45.Df 



I. INTRODUCTION 

One of the most fascinating aspects of the Anderson 
localization-delocalization transition is the occurrence of 
multifractal fluctuations of the wavefunction intensity at 
the critical point .filfil While, strictly speaking, the fluc- 
tuations are truly multifractal only at the critical point, 
where the correlation length ^ diverges, multifractal fluc- 
tuations nevertheless persist on either side of the tran- 
sition on length scales less than the correlation length.'^ 
The persistence of the fluctuations can be clearly seen in 
Fig. Ill where the wavefunction intensities for some typi- 
cal critical, metallic and localized wavefunctions are plot- 
ted. In this paper we show how to exploit this persistence 
by combining multifractal analysis with finite size scaling 
to arrive at a very powerful method for the quantitative 
analysis of the Anderson transition, or any continuous 
phase transition that exhibits multifractal fluctuations: 
Multifractal Finite Size Scaling (MFSS). 

Below we describe in detail the MFSS procedure and 
we demonstrate its potential by employing it to make 
a more comprehensive analysis of the Anderson transi- 
tion in three dimensions than given in our Letter.^ As 
we describe, the MFSS procedure permits the simultane- 
ous estimation of the usual critical parameters, such as 
the location of the critical point and the critical expo- 
nent, and the multifractal exponents. Moreover, MFSS 
offers the opportunity to examine the consistency of the 
estimates of the critical parameters against the choice of 
multifractal exponent used in the scaling analysis. 

The organization of the paper is as follows. In Section 
[TTl we describe briefly the Anderson model of a disordered 
systems and the numerical simulation of the Anderson 
transition. In Section |III[ we deflne the generalized mul- 



tifractal exponents (GMFEs) used in MFSS, and derive 
the corresponding scaling laws. The relation between 
the GMFEs and the scaling properties of the probability 
density function (PDF) of wavefunction intensities is dis- 
cussed in Section IIVI In Section |V| we demonstrate the 
necessity of avoiding correlations in the wavefunctions. 



In Sections VI and VII we present results from standard 



and multifractal FSS, including estimates for the critical 
parameters and multifractal exponents. Finally, details 
of how to account for the inevitable correlations in dif- 
ferent coarse grainings of the same simulation data, as 
well as how to check the stability of the scaling fits, are 
collected in the Appendices. 



II. THE EIGENSTATES OF THE ANDERSON 
MODEL 

We consider the three-dimensional (3D) Anderson 
Hamiltonian in site basis, 



H = 5]e. |0(^l + E l*)01' 



(1) 



ihj) 



where site i = (x, y, z) is the position of an electron in 
a cubic lattice of volume L^ — where L is measured 
in terms of the lattice constant — , (i,j) denote near- 
est neighbors, periodic boundary conditions are assumed 
and Si are random on-site energies uniformly distributed 
in the interval [-1^/2,14^/2]. The L^ x L^ Hamiltonian 
is diagonalized in the vicinity of the band centre E — 
for different degrees of disorder W, close to the criti- 
cal value Wc ~ 16.5 where the localization-delocalization 
transition occurs. The eigenstates ^ = X^i ^il*!,^'^^ nu- 
merically obtained using the Jadamilu library.EHU 




FIG. 1: (color online) Wavefunctions for the 3D Anderson model near the band center {E = 0) for a system of size L'^ — 120'^. 
From left to right, we move from the metallic phase through the transition to the insulating phase, as the degree of disorder 
increases: (a) W — 15, (b) W = 16.5, (c) W = 18. Sites contributing to 98% of the wavefunction's norm (from large to small 
values) are shown as cubes whose volume is proportional to \tpi\^- The color and opacity of the cubes is chosen according to the 
value of — log^ IV"*!^- The top plot shows the cumulative norm of the wavefunction as a function of the cut-off value considered 
for — log^ l^i|^- Vertical dashed lines mark the minimum and maximum values of \tpi\^ occurring in the wavefunction plots, 
which are (a) 9.4 x 10"* ^ IV'iP ^ 0.0035, (b) 1.4 x 10"'' sC \^^\^ sj 0.015 and (c) 1.0 x 10"" ^ \jpif < 0.18. The opacity and 
color scales are indicated below the cumulative norm plot. 



TABLE I; Average number of uncorrelated wavefunctions (A/") 
considered for each disorder W for each system size L. The 
maximum and minimum numbers of states for a given W are 
shown in brackets for each L. A total of 17 disorder values in 
the interval [15, 18] were considered. 



L (A/") (A/"max,A/"min) 



20 5138 
30 5079 
40 5168 
50 5042 
60 5027 
70 5032 
80 5028 
90 5083 
100 5024 
110 4331 
120 3103 



(5006, 5374) 
(5011, 5143) 
(5012, 5351) 
(5005, 5125) 
(5009, 5082) 
(5010, 5058) 
(5013, 5048) 
(5006, 5328) 
(5020, 5041) 
(4214, 4589) 
(3000, 3757) 



We have considered only a single eigenstate per sam- 
ple (disorder realization), namely, the eigenstate with 
energy closest to E — 0. This is costly in terms of 
computing time, but as we shall show later, absolutely 
essential to avoid the strong correlations that exist be- 
tween eigenstates of the same sample. System sizes range 
from 20^ to 120'^, and disorder values are in the inter- 
val 15 ^ W ^ 18. For each size and disorder combi- 



nation, we have taken at least 5000(3000) samples for 
L < 100(> 100), for a total of ~ 904000 wavefunctions. 
The average number of states considered for each L, W 
pair is indicated in Tablejlj In Fig.jTjwe show some wave- 
functions for L = 120 around the critical point. 



III. SCALING LAWS FOR GENERALIZED 

MULTIFRACTAL EXPONENTS AROUND THE 

CRITICAL POINT 

A. Multifractality at the critical point 

It is known that the eigenstates of the 3D Anderson 
mode l (IT]) ex hibit multifractal fluctuations at the critical 
point .ll^^liU Here, we recapitulate briefly the basics of 
multifractal analysis. 

To analyze the multifractal properties of wavefunctions 
in d dimensions of a system of size L, we coarse grain the 
wavefunction intensity on a scale I < L. The system is 
partitioned into {L/l)'^ boxes of volume l'^. A probability 



A*fc 



jGbox A: 



(2) 



is defined for each box k. It is more convenient to work, 
not directly with the box probability /i, but with a related 
random variable a, defined by 



a . 



In/i 
h^A' 



(3) 



Here, A is the ratio of the box size I to the system size L 



-I- 



(4) 



Multifractahty means that, if we count the number of 
boxes N{a) for which the value of the random variable 
a falls in a given small interval [a, a + Aa], this number 
scales with A as 



N{a) - A' 



-/(") 



(5) 



in the limit that A — ^ 0, i.e. that these boxes form a 
fractal with a fractal dimension f{a) that depends on a. 
The set of all fractal dimensions /(a) is known as the 
multifractal spectrum. 

Generalized inverse participation ratios (GIPR) or q- 
moments are obtained by summing over the boxes 



R„ 



E^l- 



(6) 



For later use, it is also convenient to define 
dR„ _ 

k 



b„ 



dq 



^^^In^fc. 



(7) 



As a consequence of multifractahty it can be shown that 
at the critical point the GIPRs obey the scaling law 



{R,)-y\ 



(8) 



in the limit that A — ?► 0. Here, the brackets denote an 
ensemble average. The mass exponents r^ depend non- 
linearly on q. They are conveniently expressed in terms 
of anomalous scaling exponents A^ , 



Tq = d{q - 1) + A, 



(9) 



The multifractal spectrum and the exponents Tq are re- 
lated via a Legendre transformation. 



aq = dTq/dq, fq = qOg - Tq, 



(10) 



which defines singularity strengths aq and a singularity 
spectrum fq. 

Multifractal exponents can also be defined from the 
scaling law corresponding to the geometric or typical av- 
erage of the GIPRs 



exp(lni?q) ^ A"^" 



(11) 



The relation between the typical and ensemble aver aged 
multifractal exponents is now well understood,'^'^ and 
the multifractal properties of the 3D Anderson transi- 
tion have been thoroughly studied using this standard 
formalism .'SliSHlIl 



B. Multifractal behavior in the vicinity of the 
critical point 

To extend multifractal analysis beyond the critical 
point we define a generalized mass exponent 



Tq{W,L,l)=\n{Rq)/\nX. 



(12) 



Here, the tilde is used to emphasize that this equation 
applies throughout the critical region not just at the crit- 
ical point. This generalized mass exponent becomes the 
usual mass exponent Tq at the critical point Wc in the 
limit A — > 0. 

Next we proceed to suggest a finite size scaling law 
for the GIPRs and from that derive a scaling law for 
these generalized mass exponents. Close to the transi- 
tion, we may suppose that the GIPRs are determined 
by the ratios of the length scales I and L to the lo- 
calization (correlation) length in the insulating (metal- 
lic) phase ^. This can be justified using renormalization 
group arguments, ^ an d is the basis for the scaling the- 
ory of localization.lSllolxhis leads to the following scaling 
law for the GIPRsP 



(i?,)(I^,i,/) = A^'7^Ji/e,VC)■ 



(13) 



At the critical point Wc the correlation length has a 
power law divergence 



^x\W -Wc 



(14) 



described by a critical exponent v. Thus, the scaling 
correctly reproduces the invariance of the GIPRs with A 
exhibited by Eq. ([s]) at the critical point. 

The scaling law for the GIPRs can be rearranged as 
follows to give a scaling law for the generalized mass ex- 
ponents. 



?q(W,L,l)=Tq 



g(g-l) 

In A 



Tq{L/tl/0. (15) 



The function Tq is related to the original TZq. The factor 
q{q — 1) has been explicitly included so that 



To = Tq = -d, Ti=Ti= 0. 

The generalized anomalous scaling exponents 



will then obey 

A,(W^,i,0 = A, 



d{q-l), 



g(g-i) 

In A 



Tq{LlUli)- 



(16) 



(17) 



(18) 



By exact analogy with Eq. ( 10 ) we may define gener- 
alized singularity strengths 



dTq/dq^ {Sq)/{{Rq)\x,\). 



(19) 



The scaling law for these quantities follows immediately 
from Eq. (flSl), 



aq{W,L,l) 



1 



Aq{L/U/0^ 



(20) 



We may also define a generalized singularity spectrum 



fg = qaq - Tq, 

obeying a corresponding scaling law, 



(21) 



fqiW, L, l)^fq+ J^^q (L / ^ I/O ■ (22) 

The scaling law for the GIPRs can be expressed in the 
entirely equivalent form 



(i^,)(W^,L,A) = A^'7^,(L/e,A) 



(23) 



This remark applies equally well to the scaling laws for 
other quantities. Writing the scaling laws in this way 
immediately suggests a standard FSS analysis by fitting 
the disorder and system size dependence at fixed A. As 
we show below this is indeed possible and works well. It 
does not, however, permit the estimation of the various 
multifractal exponents. 

A much more exciting application of the above scaling 
laws ( 15 )-( 22 ) is to fit the variation with disorder, system 



size and box size. This allows not only the estimation 
of the usual critical parameters Wc and v but also the 
simultaneous determination of a multifractal exponent 
for a particular q. Moreover, the use of different moments 
q of the wavefunctions and different averages (ensemble, 
typical) provides a test of the stability of the estimates 
for the critical parameters, as these should be average- 
and g-independent. 



C. Averaging and the A — > limit 

We emphasize that the generalized multifractal expo- 
nents (GMFEs) Tq, Aq, (Xq are equal to the corresponding 
scale invariant multifractal exponents r^, Ag,ag only at 
the critical point W = Wc in the limit A — >■ 0. 

Recalling the trivial scaling of Rq and i?i , we see that 
GMFEs Tq = — d and ti = 0, independent oiW, L, and 
I. This is equivalent to Aq = Ai = 0. 

To understand better the behavior of the GMFEs as 
functions of q, which is depicted schematically in Fig. [2] 
we discuss the expected behaviors in the thermodynamic 
limit in the metallic and insulating phases, and at the 
critical point. 

In the metallic phase, for sufficiently large system 
size or sufficiently small disorder, the states are ho- 
mogeneously extended and /ife —>■ A''. It follows that 
d{q — 1) and A^ — > for all q. 



Rq -^ A'^^'J-i) 



In the insulating phase, for sufficiently large system 
sizes or large enough disorder, the wavefunctions ap- 
proach an extremely localized state and fik -^ Sk.ko- 



Then, for q > 0, R„ ^ 1, t„ 



and Aq -> 



-d{q 



-oo 



While for g < 0, the moments diverge, t, 
and Aq — >■ — oo. We emphasize that the insulating 
limit has been confirmed by analytical calculations for 
3D exponentially localized states (with finite localization 



lengths) ^ We note that the metallic and insulating lim- 
its can be reached either as A — > or as L — > oo for any 
fixed A. 

At the critical point, we recover the multifractal expo- 
nents as A — >■ 0. These are known to obey a symmetry 
relation,^- 



A„ = Ai. 



or equivalently 



an + a 



1-9 



= 2d. 



(24) 



(25) 



This symmetry has been studied for Anderson transitions 
in different systems and dimensionality,'!^^!}!!!! qj^(\ ]-^a,s 
also been experimentally measured.'^ The corresponding 
limits for ciq can be obtained following similar reasonings 
and calculations.'22' 

in terms of the 



We may also define t^^p , AfP , ^fP , 



typical average for the moments Rq (see Eq. (11)). The 



scaling laws have the same form as for the ensemble av- 
erage exponents. For the benefit of the reader, we collect 
in Table II the definition of the exponents A^ and Sg for 
both cases. 



IV. SCALING FOR THE PDF OF 

WAVEFUNCTION INTENSITIES AROUND THE 

CRITICAL POINT 

It is also possible, and indeed, sometimes more conve- 
nient, to work directly with the probability density func- 
tion (PDF) 'P{a;W,L,l) = r{a) of a.^° At the critical 
point the distribution is multifractal 



V{a;W = Wc,LJ) oc ^/\\nX\X 

A- — ^0 



d-f(a) 



(26) 



and exhibits scale invariance, provided that A is held 
fixed. Equation (26) can be used to estimate the mul- 



tifractal spectrum directly from numerically calculated 
histograms of a values.l^ When we broaden attention to 
the critical regime, we find that, when A is sufficiently 
small, we may approximate the PDF using the general- 
ized multifractal spectrum (21 1 



7'(a;PF,L,0 oc v^h^A'^-/^"'^'^''). 



(27) 



This relation could be the basis for an alternative defi- 
nition of the generalized singularity spectrum from the 
PDF. While there will be quantitative differences with 



(21 ) we would expect the results of scaling analysis to be 



unchanged. 

The GMFEs can be obtained from the PDF. We con- 
sider So as an example. Setting q — in Eq. ( 19 ) we 
have 



ao 



A-'' In A ■ 



(28) 



TABLE II: Definitions of GMFEs for ensemble and typical averages. For a particular {W,L,l), Rq and Sq are calculated for 
each wavefunction and the ensemble average (■ • ■ ) taken over samples. Formulae for error estimation are also given, where 
a stands for the standard deviation. Note that (Sq) and (Rq) are highly correlated and their covariance must be taken into 
account in the error estimation. Definitions of the ensemble average exponents in terms of the PDF 'P{a) are also given. 



GMFE 



Error estimation 



A^^'4^-dig~) 



Definition using the PDF 



AfP 



In A 

(InRq) 

In A 
{Rq)lnX 



d{q — 1) cr^typ 



In A \R, 



(Rq) In \ 

^ 0'(lnfl,> 

In A 



1 
h^A 



A, = -dq + 



2 /C \2„2 



{R,}^ 



{R,Y 






K{Sq},{R,)) 5,= 



j^ln/„°°A«"P(a)da 

J^ a\'>°'V{a)da 
j^ \i°'V{a)da 



cr~typ — 



In A 





1. 






/- 


INSULATOR 




















N 


.\J METAL 


/'1 







... T 1 .... 


2d 




^^\ \ 


«, 


\ 


d 


\ METAL 


i \CRITICAL 





\ \^ 


? i INSULATOR 



1 



1/2 1 



1/2 

q 



FIG. 2: (color online) Schematic phase diagrams for the GMFEs in the limit A — >■ for the metallic {W < Wc), critical 
{W = Wc) and insulating (W > Wc) regimes. The metallic and insulating limits can also be reached as L — >■ cxa at fixed A. 
However, the multifractal exponents are obtained only at the critical point as A — !> 0. 



The result of averaging is the same for all boxes, so 



ao 



[a 



a'P{a)da. 



(29) 



Thus 3o corresponds to the mean value of the PDF. Ex- 
pressions for general q are given in Table III] 

The scaling of 7^(0!) with system size for fixed A in the 
vicinity of the transition is shown in Fig. l3] This figure 
shows clearly that the scaling of the distribution of wave- 
function intensities, or the distribution of a related quan- 
tity such as the local density of states (LDOS), could be 
used to characterize the Anderson transition.EHlndeed, in 
Ref.[6]the critical parameters were successfully estimated 
from a scaling analysis of the system size and disorder de- 
pendence of the maximum of the PDF of a at fixed A. 
A similar procedure might be applied to experimental 
LDOS data obtained using STM technique s,'^ ** ^ M. or 
by direct imaging of ultracold atom systems.'^^'^ 



V. REMARKS ON CORRELATIONS AMONGST 
WAVEFUNCTIONS 

In exact diagonalization studies of disordered systems, 
it is common practice to average over eigenstates located 



within a small energy window. This is because the ini- 
tialization and diagonalization of very large matrices is 
computationally demanding. Time can be saved by gen- 
erating several eigenstates for the same sample. There 
is, however, a price to be paid. Eigenstates of the same 
sample are correlated, since they are solutions of the 
Schrodinger equation with the same potential. We have 
found that these correlations distort the statistical analy- 
sis. In particular, the error estimation becomes unreliable 
and the precision of critical parameters is overestimated, 
i.e. the error bars are erroneously small. 

To quantify the correlations amongst wavcfunctions 
from the same sample, we have studied the average corre- 
lation between two eigenstates with energies E, E' , close 
to £■ = 0, such that \E - E'\ < 1. The correlation is 
defined as. 



C{q,l) 



C0v(/l9,^'«) 

a(/i9)a(^'9) ■ 



(30) 



Here, /j,' is the q-th power of the box probability Eq. (pi) 
for the eigenstate with energy E, and the prime indicates 
the same quantity calculated for the eigenstate with en- 
ergy E' . The covariance is calculated using (B2|, where 
the sum is over boxes. The covariance is normalized to 
the product of the standard deviations, where again the 



L=100 
1 = 30 
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FIG. 3; (color online) Reproduced from Ref. El Evolution 
of the wavefunction intensity distribution 'P{a;W,L,l) as a 
function of disorder W across the Anderson transition, at 
fixed A — 0.1 for two system sizes L. Each distribution was 
computed with 10'' wavefunctions. The data points (•) and 
solid lines on the bottom plane mark the trajectories of the 
maximum. For clarity, distributions are shown at W — 15, 
16.6 and 18 only. 




FIG. 4: (color online) Averaged correlation (over 1000 sam- 
ples) between eigenstates as a function of q, for L^ = lOO'' and 
varying box-size I and disorder W. Data on the left were ob- 
tained using pairs of eigenstates from the same sample. Data 
on the right, for eigenstates of two different samples. Note the 
different scales for the ordinate axis. Whenever not shown, 
errors are smaller than the symbol size. 



sum is over boxes. The correlation can thus be calculated 
for two eigenfunctions from the sample, or two eigenfunc- 
tions from a pair of samples. A further ensemble average 
over 1000 samples, or pairs of samples, as appropriate, 
was taken to arrive at Fig. |4j On the left is the corre- 
lation C{q,l) obtained using pairs of eigenstates of the 
same sample. On the right is the correlation obtained 
using eigenstates of two different samples. A maximum 
correlation (anticorrelation) corresponds to C{q,l) = 1 
(—1), whereas statistical independence implies a vanish- 
ing C{q,l). The correlation is shown as a function of 
the moment q, for different degrees of coarse-graining I 
and for three values of the disorder W = 15, 16.5, 18. 
(Note that, in the limit g — > 0, the numerator and the 
denominator in ( |30[ ) both go to zero. The finite solution 
of this indetermination is not necessarily unity, as one 
may naively expect.) Near the critical point, the corre- 
lation between eigenstates of the same sample is high, 
while for eigenstates of different samples correlation is, 
as expected, absent. We conclude that it is not safe to 
include in the MFSS analysis more than one state from a 
given sample, as these correlations render the subsequent 
statistical analysis unreliable. 

On the left of Fig. |4] we observe how the correlation 
varies with q. For large positive q, large wavefunction 
amplitudes dominate the q-th moment of the box prob- 
ability, while for large negative q small amplitudes dom- 
inate. In both cases the correlation is reduced. We 
think that this occurs because only a small number of 
boxes contribute significantly to the extreme values of 
the distribution. We have checked that the overall shape 
of the correlations in the insulating regime agrees with 
calculations for ID exponentially localized states. An- 



other feature is that for negative q correlations are ab- 
sent for box size / = 1 but are restored when the wave- 
function is coarse-grained. This may indicate that very 
small wavefunction amplitudes are affected by random 
noise, which averages away when the wavefunction is 
coarse-grained. We have verified that similar behavior 
occurs for wavefunctions calculated using other numeri- 
cal libraries, e.g. LAPACK.^ We strongly recommend 
the use of coarse-graining to evaluate negative moments, 
even for 1-D models, such as power-law random banded 
matricea^Sl and others .'231 

We emphasize that in this work — with the exception 
of the data on the left of Fig. |4] — we have used only one 
eigenfunction per sample. 



VI. SINGLE-PARAMETER SCALING AT 
FIXED A 



We study first the scaling of the GMFEs A^, 
v'yp at fixed A = 1/ lW This simplifies the scaling laws 



A*/P, 5, 



(18) and (20), which become one parameter functions. 



r,(iy,L) = g,(L/o, 



(31) 



where F^ denotes any of the above mentioned exponents. 
Since A is fixed, the multifractal exponents cannot be es- 
timated, as they merge with the zero-th order term in 
the expansion of the scaling functions Tq and Aq (see 
Eqs. (18) and (20)). Nevertheless the critical disorder 
Wc and the critical exponent v can be estimated. We 
demonstrate the consistency of the estimates of these pa- 
rameters for different q values, and their independence on 
the type of average (typical or ensemble) considered. 



A. Expansion in relevant and irrelevant scaling 
variables 



In order to fit data for the GMFEs, we follow the stan- 
dard procedure of Ref.UTjand include two kinds of correc- 
tions to scaling, (i) nonlinearities of the W dependence of 
the scaling variables, and (ii) an irrelevant scaling vari- 
able that accounts for a shift with L of the apparent crit- 
ical disorder at which the rq{W,L) curves cross. After 
expanding to first order in the irrelevant scaling variable, 
the scaling functions take the form 

g,{gL'/'',rjLy) - g'^igL^/n + vL^'QligL'/n- (32) 

Here, g and rj are the relevant and irrelevant scaling vari- 
ables, respectively. The irrelevant component is expected 
to vanish for large L,soy < 0. Both the scaling functions 
are Taylor-expanded 






G.igL' 



2_^ako0 

3=0 



'^L^/", forfc = 0, 1. 



(33) 



The scaling variables are expanded in terms oiw = {W - 
Wc) up to order irig and m,^, respectively. 



g{w) 



m=2 



b,y,W''' 



ri(w) 



rri— 1 



c„,w 



(34) 



The fitting function is characterized by the expansion or- 
ders TiQ, ni, TTig, 771,,. The total number of free parameters 
to be determined in the fit is Np — rig +ni+rng + TO^ -1-4 
(including v, y and Wc)- 

The localization (correlation) length, up to a constant 
of proportionality, is ^ = \g{w)\^" . After subtraction of 
corrections to scaling 



T,{W,L)~r^Lygl{gL^'^), 



(35) 



and the data for the GMFEs should fall on the single- 
parameter curves 



g^(±{Lio"'')- 



pcorr 

9 ~ ^q 



B. Numerical procedure at fixed A 



(36) 



When performing FSS, the aim is to identify a stable 
expansion of the scaling function that fits the numeri- 
cal data. The best fit is found by minimizing the x^ 
statistic over the parameter space. The validity of the 
fit is decided by the p- value or goodness-of-fit. We take 
p ^ 0.1 as the threshold for an acceptable fit. As a rule 
of thumb the expansion orders ng, ni, rUg, rriri are kept as 
low as possible while giving acceptable and stable fits. 
Once a stable fit has been found, the precision of the es- 
timates of the critical parameters is estimated by Monte 
Carlo simulation, i.e. by fitting a large set of synthetic 
data sets generated by adding appropriately scaled ran- 
dom normal errors to an ideal data set generated from 



the best-fit model. A detailed description of the FSS 
procedure, with some examples, is given in Appendix [A] 
We performed a detailed FSS analysis for A^ and a^ for 
both ensemble and typical averages for 13 different values 
oiq G [—1, 2] at a fixed value of A = 0.1, i.e. the box-size is 
always I = L/10. The GMFEs were obtained for system- 
sizes ranging between L — 20 and L = 120 and for 17 
values of the disorder W G [15, 18]. The average number 
of independent wavefunctions involved in the calculation 
for each L, W is indicated in Table |l] 



C. Results for A = 0.1 



The details of the fits are listed in Table [mp' In Fig. [s] 
we plot the estimates of Wc, v, and j/, along with their 
95% confidence intervals, as functions of q. 

The estimates of the critical disorder, critical exponent 
and irrelevant exponent, both for different values of the 
power q and the type of average considered, are mutually 
consistent. The values are also consistent with previous 
estimates obtained using transfer-matrix methods.'*^ *'^ ^ 
The value of the irrelevant exponent is not directly com- 
parable with previous transfer matrix studies since it is 
not clear that the dominant irrelevant correction should 
be the same for wavefunction intensity data. 

In Fig. |5] there is a clear tendency for the error bars 
to increase for g > 1. This occurs because the average 
precision of the data (see the inset in Fig. [5]) degrades 
quickly for g > 1, especially for the ensemble average. 
In general, for the same number of states, the precision 
is better when typical average is considered. Also, for 
q < 1, the precision of the data for 5^ is significantly 
better than that for A^, although this doesn't translate 
into smaller uncertainties for the critical parameters. 

In Fig. [6] we show some examples of the calculated GM- 
FEs together with the corresponding best fits. The FSS 
plots can be understood from the phase diagrams for the 
GMFEs discussed in Section |III| and depicted in Fig. [2] 
The GMFEs for W <Wciy/ > Wc) tend towards the 
metallic (insulating) limit as L — >■ oo, which is located 
above or below the corresponding critical value depend- 
ing on q. We recall that, at fixed A, the value of the 
GMFE at the crossing point Wc does not correspond to 
the scale invariant multifractal exponent, which is recov- 
ered only as A — >■ 0. 



D. The range < g sC 1/2 

For < g ^ 1/2, the expected critical value of 3g is not 
located between the values of the metallic and insulating 
limits (see Fig. [2]). This anomalous behavior makes a 
reliable FSS analysis in this interval of q very difficult. 

For q — 1/2 the values in the metallic limit and at the 
critical point are expected to be the same (= d), while 
the value in the insulating limit is expected to be zero. 
As shown in Fig. JTJJa), for q = 1/2 and W < Wc the 



TABLE III: The estimates of the critical parameters, together with 95% confidence intervals, from single-parameter FSS at 
fixed A — 0.1. Results for both ensemble (ENS) and typical (TYP) averages, for different values of q € [—1, 2], are given. Notice 
that for 5o, and also for 5i, the typical and ensemble averages coincide. The number of data is Nd and its average percentage 
precision is given in parentheses. The number of free parameters in the fit is Np, x^ is the value of the chi-squared statistic for 
the best fit, and p is the goodness of fit probability. The orders of the expansions are specified in the last column. The system 
sizes used are L € [20, 120], and the range of disorder is W £ [15, 18]. 
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FIG. 5: (color online) The estimates of the critical disorder M4, critical exponent v, and irrelevant exponent y, as functions 
of g, obtained from single-parameter FSS at fixed A = 0.1. Results for different averages and different exponents have been 
slightly offset in the q direction. Error bars are 95% confidence intervals. The corresponding values are listed in Table [Till The 
inset in the center plot shows the average data precision versus q for the different GMFEs considered. 




FIG. 6: (color online) Plots of several 
GMFEs for A = 0.1 as functions of disor- 
der at various system sizes L € [20, 120]. 
The error bars are standard deviations. 
The lines are the best fits listed in Table 
|III[ The estimated Wc are shown by ver- 
tical dashed lines and 95% confidence in- 
tervals by the shaded regions. The insets 
show the data plotted vs Lj^ with the ir- 
relevant contribution subtracted and the 
scaling function (solid line). 
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FIG. 7: (color online) Behavior of (a) qq.s and (b) 5o.45 for 
A — 0.1 as functions of disorder at various system sizes. The 
dashed lines are tentative fits that serve as guides to the eye 
only. The vertical dotted lines indicate the position of Wc 
according to these tentative fits. The horizontal dotted line 
in (a) indicates the metallic limit ao.s ~ 3. The insets display 
the corrected data and the scaling functions according to the 
tentative fits. 



data are almost independent of both L and W. While, 
for W > Wc, the dependence on L is very strong, as the 
data tend to the insulating limit. This behavior cannot 
be reliably described by a power series expansion. 

For < g < 1/2 the value of 5g at the critical point 
is now expected to be larger than the values in both the 
metallic and insulating limits. Data for q — 0.45 is shown 
in Fig. JTJ^b). Very close to the critical point, a standard 
scaling behavior with opposite L dependence at each side 
of the transition is visible. However, for larger disorder, 
the L dependence is again reversed so as to approach 
the value in the insulating limit. We wish to emphasize 
that this is not a numerical artefact. Rather, analytical 
calculations for 3D exponentially localized states confirm 
the behavior. 



VII. MULTIFRACTAL ESS 



matrix and the x^-minimization procedure for correlated 
data in Appendices [B| and [C| 



B. Expansion of scaling functions 

In MFSS, the scaling functions are functions of two in- 
dependent variables L/^ and l/^. While these variables 
are independent because the system size and the box size 
vary independently, they involve the same scaling vari- 
able and renormalize in the same way. In addition, we 
need to allow for non-linear dependence in W and for ir- 
relevant scaling variables. We have found that the most 
important irrelevant contribution is due to the box-size 
I. Therefore, we use the expansion 



A,(eil/^g^l/^,7/«)- 



A„ 



Hi/L) t. 



^ir,iyfT'{QL'/-',gi'^n, 



(37) 



and similarly for ciq. Here, g and rj are the relevant and 
irrelevant scaling variables, with 1/v and y < the cor- 
responding exponents. Note that we expand to second 
order in the irrelevant variable. We find that this is nec- 
essary to fit the data reliably. The functions 71*^ are ex- 
panded, 



T'igL'^-'^gl'^n 



i=0 i=0 



for fc = 0, 1, 2, as are the fields as described in Eq. (34 1 



The expansion of the scaling function is then character- 
ized by the indices n'^]^,n^,n\,n],n\,n'f,rng,rnj^. The 
number of free parameters is given by 



We now consider the scaling behavior of the GMFEs 
as function not only of disorder W and system size L 
but in addition the box-size I using scaling laws ( 18 ) 



and (20). This permits the simultaneous estimation of 



the multifractal exponents and the critical parameters z/, 
Wc, y- This is a major advance over traditional multi- 
fractal analysis, where the position of the critical point 
must be estimated in a separate FSS analysis before the 
multifractal analysis. 



Fitting of correlated data 



2 
Np = ^(n 

fe=0 



l)(n 



After subtraction of irrelevant corrections we have 



(39) 



A-"- = ^,+T^{±{L/£,f'\±{llO"niHllL) (40) 

and the numerical data should fall on a common surface 
rather than a common curve as in standard FSS. 



When evaluated at fixed A, Eq. (37) leads to the FSS 



expansions given in Sec |VI| Since / ~ XL, when perform- 
ing FSS at fixed A the irrelevant correction is determined 
by the system size. 



Different coarse-grainings I for the same disorder W 
and system size L use the same set of eigenstates, which 
induces correlations among the estimates of GMFEs for 
different / and the same W and L. These correlations 
must be properly taken into account in MFSS. To do so 
we generalize the definition of x^ in the numerical mini- 
mization by including the full covariance matrix for the 
GMFEs. We describe the calculation of the covariance 



C. Results 

We have carried out the MFSS analysis on the ensem- 
ble averaged GMFEs A, for different g G [—1,2], and ctq 
for q — 0, 1. The estimates of the critical parameters and 
the multifractal exponents, together with full details of 



the fits, are included in Table IV 
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We have considered different ranges of data for dif- 
ferent q trying to maximize the number of points that 
we could fit in a stable manner, while keeping the value 
of A ^ 0.1. The minimum values of A occurring in our 
data sets are Amin = 0.017 {Imin — 2) for q ^ 1/2 and 
Amin = 0.008 (/mill = 1) for q > 1/2. Bccausc of the need 
for coarse-graining for negative q, as discussed in Section 
Til we exclude alH = 1 data for q ^ 1/2. 

The best-fit estimates for Wc, y and t/ as functions of 
q are shown in Fig. [8J The consistency of the estimates 
of the critical parameters for different q is remarkable. 
Compared with the FSS at fixed A of Section [VH the es- 
timate of the irrelevant exponent y is more stable. As q 
increases, the uncertainty for the estimates of Wc, v and 
y grow. This is partly due to the loss of precision in the 
data for high q (see the inset in Fig. |8|, but this does 
not fully explain the large difference in the uncertainty 
between sets with similar precisions, for example q = —\ 
and q = 1.5. We believe that the difference is caused by 
the amplitude of the irrelevant component in the data, 
which is larger for q < 1/2. (We have confirmed this by 
examining the q dependence of the cocfBcicnt aioo in the 
expansion (37).) If the amplitude of the irrelevant shift 



in the data is small, the estimation of the irrelevant expo- 
nent becomes very ambiguous. This in turn leads to an 
increase in the uncertainties of all the other parameters. 
The best precision is achieved for Sq (0.05%). From the 
MFSS analysis of this GMFE we find 



and 



Wc = 16.530(16.524,16.536) 



V = 1.590 (1.579, 1.602) 



(41) 



(42) 



where the error limits correspond to 95% confidence in- 
tervals. 



D. Scaling surfaces 

In Fig. [9] we show the best fits and the correspond- 
ing scaling surfaces for A_o.75 and cxq. The upper plots 
display the GMFEs and cross-sections of the global fit 
displayed at the different A values occurring in the data 
set. The bottom plots show the scaling functions in terms 
of the two variables L/^ and A. (The visualization of the 
scaling functions is improved when A is chosen instead 
of l/£, as the second variable.) The scale invariant mul- 
tifractal exponents correspond to the asymptotic value 
at the critical point as A — > 0. This is highlighted in 
the insets of Fig. [9lc,d), where the behavior of the scal- 
ing function at criticality — when the sheets of extended 
and localized phases meet — is shown versus log(A). 

The transition can also be visualized as shown in 
Fig. [To] where we display a cross-section at fixed box- 
size of the data and the fit for A_o.75 and 5i. This 
provides an alternative way to monitor the Anderson 



transition.ESHlZl \\/q see the flow with increasing L (de- 
creasing A) towards the metallic or insulating limits de- 
pending on the disorder W, and at the critical point Wc 
the convergence to the multifractal exponent. At the 
critical point, the scaling law (37 1 reduces to 



1 2 



(43) 



k=0 



From this we find that this convergence differs for differ- 
ent I. This is illustrated in the insets of Fig. [TO} 



E. The multifractal spectrum 

Previous work by some of the authors, using standard 
multifractal analysis at W — 16.5 ~ Wc, suggested that 
the symmetry relatio n (24) for Ag holds for the 3D An- 
derson model.'SI^SI^Sl Our improved analysis given here 
confirms this at Wc- As shown in Fig. [TT] the expo- 
nents Aq satisfy the symmetry relation for the range of 
q values considered. The symmetry relation ( |25[ ) for a^, 
is also satisfied by our estimates ao = 4.048(4.045, 4.050) 
and ai — 1.958(1.953,1.963). A more careful analysis 
in terms of the reduced exponents Aq/q{l — q) reveals a 
slight violation of the symmetry beyond q = 1.75. We 
suspect that this reflects minor limitations in our A — > 
extrapolation, rather than a genuine violation of the sym- 
metry relation. (Particulary for q < where the range 
of A is more restricted because of the absence of data 
for 1 = 1.) We note that the validity of the symmetry 
relation for the whole g-range depends on the absence 
of termination points for the ensemble average multifrac- 
tal spectrum,'^' which is still an open question for the 3D 
Anderson transition.ESl 

Our data in Fig. [Tl| (inset) clearly show that the mul- 
tifractal spectrum is not parabolic, i.e. that the reduced 
exponents depend on q. This is in agreement with pre- 
vious results at the 3D Anderson transition, where non- 

parabolicity has been directly observed in / (g),-^ and in 

^[3oiT 



the non-Gaussian nature of the PDF, V{a)^^^^A non- 
parabolic multifractal spectrum implies that the distri- 
bution of wavefunction intensities at the transition, and 
consequently also of the LDOS, is not a log-normal dis- 
tribution, in contrast to recent claims.'^ 



VIII. CONCLUSIONS 

We have shown how to exploit the persistence of 
multifractal fluctuations away from the critical point 
on scales below the correlation length to perform a 
multifractal finite size scaling analysis. We demon- 
strated the potential of this approach by applying it 
to the Anderson localization-delocalization transition in 
the 3D Anderson model. We validated our proposed 
scaling laws for the generalized multifractal exponents 
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TABLE IV: The estimates of the critical parameters and multifractal exponents together with 95% confidence intervals, from 
MFSS of Ag for q £ [—1, 2] and a.q for g = 0, 1, under ensemble average. The number of data used is Nd (average percentage 
precision in parentheses), the number of free parameters in the fit is Np, x^ is the value of the chi-squared statistic for the best 
fit, and p is the goodness of fit probability. The last column specifies the orders of the expansion: n1,n'^,n]^,nl , n\,nf, nig, nir,. 
The system sizes considered are L £ [20, 120], the range of disorder is W & [15, 18], minimum box-size Imin = 2 (Amin ~ 0.017) 
for q ^ 1/2 and /min ~ 1 (Amin = 0.008) for q > 1/2. The maximum values considered for A change from Anmx = 0.063 to 
Amax = 0.1 for different q. 
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FIG. 8: (color online) The estimates of the critical parameters v, Wc and y, as functions of q, obtained fr om M FSS for A^ and 
3q (only q = 0,1). Error bars are 95% confidence intervals. The corresponding values are listed in Table IV The inset in the 
centre plot shows the average data precision versus q for the data set used. A density plot of the histograms obtained from the 
Monte Carlo simulations used to determine the uncertainty of the estimates is shown for each q. The color scale on top of each 
graph is for the density plot. The histograms are normalized so that their maximum value is unity. 



and we estimated the critical disorder at the band cen- 
ter Wc = 16.530(16.524, 16.536) and the critical expo- 
nent describing the divergence of the localisation length 
ly — 1.590(1.579,1.602). A remarkable consistency of 
these estimated critical parameters for different averages 
(ensemble and typical) and different generalized multi- 
fractal exponents (different q) was seen. We also found 
that the multifractal spectrum exhibits the predicted 
symmetry and confirmed its non-parabolic nature. 

Recently, the value of v for the Anderson transition has 
been estimated for a variety of systems, including disor- 



dered phonons^, electrons in non-conventional lattices^ 
or with topological disorder,'S3 a disordered Lenard- Jones 
fluid,'^ and cold-atom quantum kicked-rotors, '"'''' '''^ find- 
ing agreement with our results. The universality of the 
critical exponent for the Anderson transition is, there- 
fore, well established and a value for i/ ~ 1.6 is widely 
accepted by the community. 

In our approach, scaling versus system size yields the 
critical exponent, while scaling versus box size yields the 
multifractal exponents. Thus, these exponents do not ap- 
pear to be intrinsically related. Nevertheless, such a con- 
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FIG. 9: (color online) MFSS of A_o.75 (a,c) and 5o (b,d). (Upper plots) GMFEs (•) as functions of disorder W, for different 
A = l/L. The solid lines are cross-sections at fixed A of the best fit, plotted for different L. Note that all points are fitted 
simultaneously. Alternating colors have been used for better visualization. (Lower plots) GMFEs with irrelevant contribution 
subtracted (•, o) and the scaling surfaces (symbol o highlights the maximum value of A). The insets are the scaling functions 
at the critical point, highlighted also in the right face of the main plot. The arrows indicate the multifractal exponents given 
by the extrapolation A — >■ 0. The shaded regions indicate the range of A accessed in our simulations. 



nection has been suggested by Kramer^ and Janssenjl^^l 
who discussed a lower bound for v in terms of the multi- 
fractal exponents, namely 



v > 2/t2. 



(44) 



From our estimate of A2, the confidence interval for 2/t2 
is (1.571, 1.656). Remarkably, all our estimates for 1/ fall 
within this interval. Although this bound has not been 
widely discussed in the literature, we find the agreement 
intriguing. 

Multifractal fluctuations have been observed in the 
local density of states near the metal-insulator transi- 
tion in semiconductor s'^ '^^J and in the intensity fluctu- 
ations of ultrasound ne ar th e Anderson transition in a 
random elastic network.l^ES We believe that multifrac- 
tal finite size scaling as demonstrated here has potential 
in the quantitative analysis of this sort of experimental 

data.'S2n36l55i56i59] 
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Appendix A: Fitting procedure 
1. Table of fits. 

The first step is to generate a table of fits by vary- 
ing the expansion indices of the model over a reasonable 
range of values. Since the number of combinations of 
the indices increases exponentially, we impose sensible 
restrictions. Let us consider, for example, the single- 



parameter FSS of Section VI where the expansion ( 32 ) 
is characterized by the indices no, ni, nig, 



m^. We ex- 
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FIG. 11: (color online) Multifractal exponents A, obtained 
from MFSS. The numerical values are listed in Table IIVI 
Error bars denote 95% confidence intervals, and are con- 
tained within symbol size whenever not shown. Note that 
Ao = Ai = by definition. The inset shows the reduced 
multifractal exponents Aq/q{l — q). The horizontal dotted 
line corresponds to Wegner's parabolic approximation.'^ 



fit. We use p ^ 0.1 as a criterion for an acceptable fit. 

After excluding unacceptable fits we order the table 
according to tlie number of free parameters. 



FIG. 10: (color online) Plot of (a) A_o.75 and (b) 5i for fixed 
box-size I as a function of A and increasing disorder values 
W e [15, 18] from top to bottom. The lines are cross-sections 
at (a) Z = 2 and (h) I — 1 of the global multifractal FSS fit. 
The thick lines (red ) correspond to the GMFE at the critical 
point Wc [Eq. ( |43| l]. The insets show the latter function for 
different values of I. The dashed line in the inset displays 
the GMFE at criticality with vanishing irrelevant corrections 
(effectively I = oo). The circle on the inset y-axis marks the 
asymptotic (A — >■ 0) multifractal exponent. 



pect the irrelevant components to be less important than 
tlic relevant part of the scaling function. Therefore we 
impose the restriction ni ^ tiq and m,, ^ rUg. For the 
single-parameter FSS we generate tables for each q con- 
taining usually a few hundred index-combinations. For 
the MFSS of Section |VII| we generate initial tables of 
~ 2000 index-combinations for each q and for different 
sets of data. 

We decide whether or not a fit is acceptable based on 
the goodness-of-fit or p-value. This is the probability 
that for k degrees of freedom we would observe a value 
of x^ larger than the one obtained for the best fit. To 
calculate this probability we use the approximation 



P 



r(fc/2,xV2) 
r(fc/2) ■ 



(Al) 



Here, T{x) is the Euler gamma function, T{a,x) is the 
upper incomplete gamma function, and k denotes the 
number of degrees of freedom in the model, i.e. the num- 
ber of data minus the number of free parameters in the 



2. Stability of fits. 

A good fit must not only give an acceptable p value, 
it must also be stable. By stable we mean that the esti- 
mates of the fit parameters, and specially of the critical 
parameters, must not change significantly when the ex- 
pansion orders are increased. In order to check the stabil- 
ity of a particular expansion, for example hq, tii, mg, rriri 
in fixed- A FSS, we consider all expansions where each in- 
dex is separately increased by 1, and also the case where 
all indices are increased by 1 at the same time. If the 
critical parameters from all these fits lie within the un- 
certainty interval obtained for the estimates of the initial 
expansion, then we regard that fit as stable. The search 
for a stable fit proceeds by considering acceptable fits 
from the initial table in increasing order with the num- 
ber of free parameters. As a general rule we try to find 
the simplest stable fit. We give examples of this proce- 
dure in Tables |V] and I VII While the criterion for a stable 
fit can certainly be debated we think that our criterion 
is sensible and helps avoid the ambiguity that may arise 
from an arbitrary choice of a particular fit with an ac- 
ceptable p value. 



3. Uncertainties for fitting parameters. 

To determine the uncertainties (confidence intervals) 
of the estimates of the fit parameters we use a Monte 
Carlo method. 
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TABLE V: Table of fits for FSS of Al^Pg at fixed A = 0.1. 
Our clioice of fit is fiigfiliglited in bold. In subsequent rows the 
expansion indices are progressively increased. Dashes indicate 
the expansions used to check if our chosen fit is stable. The 
stability of our choice is shown in Fig. |12| 
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TABLE VI: Table of fits for MESS of qq. Our choice of fit 
is highlighted in bold. The stability of our choice is shown in 
Fig.fll 
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FIG. 12: (color online) Histograms of the critical parameters 
for expansion (312 0) for FSS at fixed A = 0.1 of A^Z.^'g. 
The histograms are obtained from fits of 2000 synthetic data 
sets. The j/-axis indicates the number of events. The vertical 
dashed lines denote the 95% confidence interval. The circles 
on the a;-axis mark the position of the values obtained from 
the stability analysis shown in Table [V] The color scale is the 
same as that used for the histogram density plots in Fig. [8] 



feet data set with the same standard deviations as the 
simulation data to generate the synthetic data sets. In 
Fig. [12] we show the resulting histograms for the fit of 
A*Z^5 affixed A = 0.1. 

For MFSS the simulation data are correlated and these 
correlations must be taken into account when generating 
the synthetic data sets. In this case, the appropriate dis- 
tribution is a multivariate normal distribution centered 
on the perfect data set and with the covariance matrix 
estimated from the simulation data as described in Ap- 
pendix [Bl A simple way to generate the necessary corre- 
lated data is to use the Cholesky decomposition of the co- 
variance matrix to define an auxiliary set of independent 
normally distributed random variables with unit variance 
(see Appendix [C]) . In Fig. 13 we show the histograms for 
critical parameters obtained from the MFSS of So- 



Given a candidate expansion for a stable fit, we gen- 
erate a perfect data set from that fit. We then generate 
at least 1000 synthetic data sets by sampling randomly 
from an appropriate (see below) multivariate normal dis- 
tribution centered on the perfect data set. We then fit 
these data sets, using the same expansion of the scaling 
functions, and build histograms of the estimates for the 
critical parameters. From these histograms, the corre- 
sponding, possibly asymmetric, 95% confidence interval 
is obtained by removing 2.5% of the events at each end 
of the distribution. 

For single-parameter FSS at fixed A, the simulation 
data are uncorrelated and it is sufficient to use a product 
of independent normal distributions centered on the per- 



Appendix B: Calculation of the covariance matrix 
for correlated data 

In this section only we distinguish between the expec- 
tation value or population mean indicated by angular 
brackets, and the sample mean indicated by an overline. 
To avoid a cumbersome notation we do not make this 
distinction elsewhere in the paper. In practical calcula- 
tions, the population mean is always estimated using the 
sample mean. 

The correlation between two random variables X and 
Y can be characterized in terms of the covariance, defined 
as 



coy{X,Y) = {XY)-{X){Y). 



(Bl) 



16 



100 




FIG. 13: (color online) Histograms of the critical parameters 
obtained for expansion (32211120) for MFSS of 5o. The 
histograms are obtained from fits of 1000 synthetic data sets. 
The y-axis shows the number of events. The vertical dashed 
lines indicate the 95% confidence interval. The circles on the 
X-axis mark the position of the values obtained from the sta- 
bility analysis shown in Table [Vl] The color scale is the same 
as that used for the histogram density shown in Fig. [S] 



Note that cov(X, X) = var(X) = aj^, gives the 
variance for the variable X. From n samples 
{Xi,Yi), . . . , {Xn, Yn), of the pair of random variables, 
the covariance can be estimated as!^ 



cov 



i=l 

Here, the overline indicates the sample mean 

1 " 

71 ^ ^ 



(B2) 



(B3) 



and similarly for Y. We also need the covariance of the 



sample means. This follows from (Bl| 



cov(X,r) = -cov(X,Y). 
n 



(B4) 



Furthermore, we need to consider random variables 
Y^^' , . . . , Y^'^' that are in turn functions of random vari- 
ables X^^), ..., X^ such that 



w = r«(x«,...,xM), 



Y^^i = 



(B5) 



where i — 1, . . . , s. Assuming that it is reasonable to lin- 
earize the functions in the range of interest, the covari- 
ance matrix C of the random variables X'^^\ . . . , X^'^\ 



C,j =cov(x«,X(^) 



(B6) 



is related to the covariance matrix C" of the variables 

y(i),...,y('^) by 



C" = J'^CJ. 



(B7) 



Here, J is the matrix of derivatives 






(B8) 



evaluated at (X(i)), . . . , (X^). 

We now consider the application of the formulae above 
to the determination of the uncertainties of the GMFE 
Ag. Correlations only arise between values of this expo- 
nent calculated for different coarse-grainings of the same 
set of samples. Therefore, the covariance matrix for this 
quantity is block diagonal in disorder W and system size 
L. Consider two different coarse-grainings It and /j, and 



(i) 



:u) 



corresponding estimates, A, and Ag , of the values of 
this GMFE for these coarse-grainings. Application of the 
formulae above is then straightforward because Ag is a 
function of Rq only (see Table HU) . The result is 



cov(a«,A(^)) = 



I Kg jitg I 



(4'^) {R\^'^)\n\<^^)\x,m 



(B9) 



Application of the formulae above to the determination 
of the uncertainties of the GMFE 3g is more complicated 
because it is a function of Rq and S'g. The result is 



(5«,5^)) 



(i)\/c.O)\ 



(Sg )(6g' 



W\2/^UA2 



{R', 



I Kg , Kg 1 



0),.cov(i?«,4^') 



(4^ Vov (5« , r'^P) (4^ Vov (4'^ , i?«) 



{R^){R^_ 
cov(4^\4')) 



(4^V(4''^) 



(i?«)(4^"') 



/(inAWlnA^^)). (BIO) 
For g = 0, 1 the formulae are considerably simpler, 



cov 



.,,.,,„. [A"'AOI]^ov(5'-',5i") 



(5«,5(- 



0) 



cov 



InAWlnA(j) 



In AW In AO) 



(B12) 



Appendix C: x^ minimization for correlated data 

Let Yi be the result of the i-th simulation performed 
for simulation parameters Xi. The uncertainties and cor- 
relations amongst the simulation data are described by 
a covariance matrix C = {cov(yi, y,)}. We fit a model 
Fi = F(Xi, {a}) to the simulation data by varying the 
parameters {a} of the model so as to minimize the chi- 
squared statistic 



X 



J2iY.-F.){C-'),AY,~F,). (CI) 
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In practice it is convenient to perform a Cholesky fac- 
torization, C~^ = R"^ R^ where R is an upper triangular 
matrix, so that x^ can be expressed as the square norm 
of a vector x^ = \\R{Y - F)||2. 



In MFSS, Yi are the values of the GMFE under con- 
sideration and Xi = {Wi,Li,li\. The covariance matrix 
is calculated as described in Appendix [B] 
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